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qv ■ Abstract. Lyapunov exponents are indicators for the chaotic properties of a classical 

CN| . dynamical system. They are most naturally defined in terms of the time evolution 

of a set of so-called covariant vectors, co-moving with the linearized flow in tangent 
Q ■ space. Taking a simple spring pendulum and the Henon-Heiles system as examples, we 
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demonstrate the consequences of symplectic symmetry and of time-reversal invariance 
for such vectors, and study the transformation between different parameterizations of 
the flow. 
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1. Introduction 

The stability of the phase-space trajectory of a dynamical system is determined by 
the so-called Lyapunov exponents, which are the time-averaged rate constants of a 
set of perturbation vectors in tangent space, which grow or shrink exponentially with 
time. Various such sets have been introduced in connection with the algorithms 
for the computation of the Lyapunov exponents. The features of these sets and of 
their associated exponents are well known to mathematicians and theoretical physicists 
[H El El H] . In view of an increasing number of applications to ever more sophisticated 
physical systems, it seems worthwhile to become familiar with the properties of these 
tangent-space objects for some simple models. This is the aim of the present paper. 

The most familiar set of perturbation vectors spanning the tangent space is the 
set of ort honor mal Gram Schmidt (GS) vectors {^g £ }, £ € {1, . . . ,D}, commonly 
also referred to as backward Lyapunov vectors (since they depend on the history in the 
past). D is the dimension of phase space. For time-continuous systems they are most 
elegantly obtained as the forward-in-time solution (indicated here by an index (+)) of 
the linearized motion equations, augmented by a set of constraints, which continuously 
enforce the orthogonality and the norm conservation of these vectors. The latter gives 
rise to the GS-Lyapunov exponents along the way |6]. Instead of constraints, the 
standard algorithms for the computation of Lyapunov exponents [7J El E] employ a 
periodic Gram-Schmidt re-orthogonalization scheme, which may also be easily adapted 
for many- dimensional systems involving time-discontinuous maps such as the dynamics 
of hard spheres [10J. 

A second and less familiar set of perturbation vectors are the covariant vectors 
{( + )i/}, £ g |x ; ... ; _D}, which are of unit length but generally not orthogonal to 
each other. They evolve in tangent space according to the non-constrained linearized 
motion equations. Still required is a periodic re-normalization, which also generates 
the corresponding Lyapunov exponents. These vectors constitute a practical realization 
of the Oseledec splitting [H El E] of the tangent space into a hierarchy of stable and 
unstable subspaces at any phase-space point visited by the trajectory. 

In addition to these sets connected with the tangent flow forward in time, there 
exist corresponding sets, if the dynamics is followed backward in time. They will be 
distinguished by the index (-). Their application gives rise to the orthonormal forward 
Gram-Schmidt Lyapunov vectors {^g 1 }, which are conventionally called "forward" 
since they depend on their history in the future. In general, they are not simply related to 
the set of backward Gram-Schmidt vectors. Similarly, there exists a set of time-reversed 
covariant vectors {^v e }, which, however, for time-reversal invariant dynamics agrees 
with its time-forward counterpart up to a simple reversal of the indices, £ — > D + 1 — £. 

While the Gram-Schmidt vectors have been central to any algorithm since the 
pioneering days of the numerical stability analysis of dynamical systems [71 El E] more 
than thirty years ago, the covariant vectors proved rather elusive due to the Lyapunov 
instability of the computational process itself. Practical schemes for the computation 
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of covariant vectors have been developed only recently [TTJ [27]. A few studies of their 
properties for various systems have appeared since then [U [131 HSl EH EE3 HEl UHl I2D] ■ 
In the following sections we shall study these properties for two simple symplectic 
systems, the chaotic spring pendulum and the Henon-Heiles system. To establish our 
notation, we first summarize the most important relations below. 

2. Definitions and notation 

If r(t) denotes the state of a dynamical system of dimension D, its evolution equation 
and formal solution are given by 

f = F(r), r(t) = 4 (r(o)) (i) 

where F is a (generally nonlinear) vector- valued function of dimension D, and the map 
(f) 1 defines the phase flow. The linearized evolution equation and the formal solution for 
an arbitrary perturbation vector ST(t) in tangent space become 

5T = J(T)ST, ST(t) = D0*|r(o) ^(0) (2) 

where J(T) = dF/dT, and where D0*|r(o) denotes the flow in tangent space. As 
already mentioned, the covariant vectors evolve - co-rotate in particular - according to 
the unconstrained tangent-space flow, 

Mr(t)) -| W | r(0) ^(r ( o))ir (3) 

Of course, the algorithm has to ascertain that the initial vector i/(T(0)) is already 
properly oriented and normalized in order to qualify as covariant. The stretching factor 
in the denominator of Equ. ([3]) provides the general definition for the (global) Lyapunov 
exponents: 

&\t= lim 1 In ||£>0*|r(o)« £ (r(O)) ||, ie{l,--,D} (4) 

t— ^±oo t 

Under very mild conditions, the multiplicative ergodic theorem of Oseledec [TJ [3] asserts 
that the global Lyapunov exponents are dynamical invariants and do not depend on the 
norm and, hence, the particular parameterization in phase space. 
Two other limits of Eq. ()3j) are of special interest: 

• r large but finite: 

A£ ov (r) = lim ^ln||L>0Vo)t/(r(O)) N r large but finite. (5) 

These objects depend on r and are referred to as finite-time Lyapunov exponents 
(FTLE). For infinite time they converge to the global exponents A^. For finite 
t they are no dynamical invariants and depend, for example, on the coordinate 
system in use. However, it was recently shown by Kuptsov and Politi [20] that the 
fluctuations of these quantities, i.e. the linear growth rate of the covariances of the 
logarithmic expansion factors tA^ ov (t), 

D™\t) = lim [(Ar(r)Ar(r)) - Ut'] T, (6) 

T<00 
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is a dynamical invariant for large-enough r. Here, (...) denotes an average over 
(infinitely) many uncorrelated realizations along a trajectory. Below we shall 
demonstrate this property for the chaotic pendulum. 

• r — > 0: The limit 

A, cov (r(O)) = lim^ln||D0^| r(o) ^(r(O)) ||, £E{l,...,D} (7) 

provides a definition of the so-called local Lyapunov exponents (LLE). They are 
point functions in phase space. Their time average over a long trajectory also 
converges to the global exponents Xi. If the covariant vector v i (T) is known at a 
point r the corresponding covariant LLE follows from 

Af v (r)=[t/(r)f^(r)t/(r), (8) 

where T means transposition and J is the Jacobian of Eq. ([2]). This nicely 
underlines the local nature of the LLEs, which depend implicitly on time. 

So far all definitions of Lyapunov exponents are in terms of covariant Lyapunov 
vectors. There are analogous definitions for the orthonormal Gram-Schmidt vectors 
both forward and backward in time, which are summarized in Ref. [16] and are not 
repeated here. The indices GS and cov will be used in the following to distinguish 
between the quantities. The GS-FTLEs also converge to the global exponents with 
t —7- oo, as do the GS-LLEs when time-averaged along a trajectory. It is interesting to 
note that the algorithm with continuously constrained orthonormality mentioned above 
El HE] provides an expression for the GS-LLEs, 

Ap s (r)=^(r)] T .7(r)g^r), (9) 

which is analogous to Eq. ([8]) for the covariant case. 

The classical algorithms involving GS re-orthonormalization keep track of the 
volume changes of <i-dimensional volume elements in phase space (d < D). For 
symplectic systems, for which the phase volume is invariant, the GS-LLEs show the 
symplectic local pairing symmetry J2TJ [161 02] , 

«Ap S (i)= -<+> AgJL(t) (10) 

H A? S W = - H AgL(t), (11) 

where ^ and ^ indicate whether the trajectory is followed forward or backward in 
time. Since the angle information is discarded by the GS-process, these exponents do 
not show time-reversal symmetry, 

HAp S (t) ^ -«AgS_^). (12) 

If the system is not symplectic, also Eqs. (TTUl) and (fill cease to exist. On the other 
hand, the orientation of the covariant vectors is not affected by any process such as re- 
orthogonalization, and the covariant LLEs display the time reversal symmetry [21 IT6] : 

(-)A, cov (r(t)) = - (+) A£° + V,(r(t)) ; e = 1, • • - , d. (13) 
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Similarly, the covariant vectors obey 

(~V(r(t)) = ±w/ +H (r(t)). (14) 

These relations apply whether or not the system is symplectic, as long as it is time 
reversible. 

Another attractive property of the covariant vectors derives from the fact that 
they are spanning sets for the stable and unstable Oseledec subspaces associated 
with the respective Lyapunov exponents (one-dimensional without degeneracy, and 
m-dimensional for multiplicity m), and are defined without reference to a particular 
parameterization. This does not apply to the GS-exponents which suffer from the 
additional constraint of orthogonality. It has been shown recently by Yang and Radons 
that covariant vectors may be easily transformed from one coordinate system to another, 
and the same is true for the covariant finite-time and local Lyapunov exponents [J] We 
shall demonstrate this property for the chaotic pendulum below by transforming from 
the Cartesian representation to a polar coordinate system. 

In the following section we consider the chaotic spring pendulum, and in Section 
3 we discuss the Henon-Heiles system. We close with a short survey of the merits and 
disadvantages of the covariant analysis. 

3. The spring pendulum 

We consider the planar chaotic motion of the spring pendulum [221 123] as a simple 
example. It consists of a mass m attached to a harmonic spring with spring constant 
k and rest length R, which is exposed to a homogeneous gravitational force of strength 
mg in the negative y direction. The spring is attached to a pivot at the origin. 
The Hamiltonian in Cartesian coordinates (x,y) and with the conjugate momenta 
(Px,Py) = m(x, y) is given by 



^2 i m 2 



Px+Py , k 



2 



H c = ^-^+ , -^x 2 + y 2 -R) +mgy, (15) 

from which follow the equations of motion for the reference trajectory in phase space, 
x = p x /m 

Px = k[(R/r)-l}x [ ) 

Py = k[(R/r) - l]y - mg, 
and for the perturbation vectors in tangent space. 

Sx = 5p x /m 

6y = 5p y /m 

5p x = k[(R/r) -\\Sx- (kRx/r 3 )(x5x + y5y) l ' 

5p y = k[(R/r) — l]Sy — (kRy/r 3 )(x5x + ySy), 



Here, r = \Jx 2 + y 2 . 
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Introducing polar coordinates (r, 0) through 

x = r sin 0, y = r cos 0, (18) 

the Hamiltonian becomes 

where p r = mf and p<p = mr 2 <j) are the conjugate momenta. The equations of motion in 
phase space now read 

r = p r /m 

<$> = P<p/(mr 2 ) 

Pr = p'i/imr 3 ) — k{r — R) — rag cos 

Ptj, = mgr sin 0, 

from which the linearized motion equations for the perturbation vectors readily follow, 

Sr = Sp r /m 

5(f) = 5p^/{mr 2 ) -2p (f> 5r/{mr i ) 
Sp r = 2p ( f > 8p<j > /(mr 3 ) — Zp\brj \mr A ) — kbr + mg sin r^ 



mg sin <pbr + mgr cos ( 

Let us consider the state point T c = (x,y,p x ,p y ), which in the polar coordinate 

system becomes T P = (r, 0,p r ,p0). This, according to Eq. (fT8|) . is accomplished by the 

transformation [4\ 

\ T 
xp x + yp y 



V(T C )= fv/^T^tan- 1 f€\ 



\/x 2 + y 1 



(p x y - xp y ) . (22) 



Any Cartesian covariant vector v c is transformed to the polar representation according 
to 

bT l P = M v £ c , (23) 

where M. = dV/dTc is the Jacobian of the transformation ( )22l) . A final normalization, 



•, = ±^ = ±1^, (24) 

yields the covariant vector in the polar representation. Although formulated in terms 
of a particular coordinate transformation, the expressions ( 1231) and (12^1) are completely 
general [1]. 

For our numerical work we use reduced units for which the mass m, the rest length 
R, and the gravitational acceleration g are unity. The spring constant k is set to two. For 
the initial condition of the reference trajectory we take (x, y,p x ,Py) = (0.00001, 1, 0, 0), 
which determines the energy. Care must be taken that the reference trajectory in the 
Cartesian and polar representations coincide. Therefore, only Eqs. ffl6|) are integrated, 
and the solutions of the polar equations ( 1201) are obtained by the Cartesian-to-polar 
transformation (1221) . The computation of the covariant vectors and their associated 
local Lyapunov exponents has been outlined in Refs. (TTJ [151 [IB], to which we refer for 
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details. The integration is carried out with a 4th-order Runge-Kutta algorithm with a 
time step of 0.002. The global Lyapunov spectrum is found to be {0.0565, 0, 0, —0.0565}. 
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Figure 1. (Color online) Projection of a short spring-pendulum trajectory onto the 
(a;,2/,p2;)-subspace. The blue and green lines represent projections of the covariant 
vectors v 1 and v 4 , respectively, and give an impression of the respective unstable and 
stable manifolds along the trajectory. The time interval for this trajectory segment is 
16 time units. 



By the smooth (red) line in Fig. [T]we show a projection of a short trajectory into 
the three-dimensional subspace spanned by x,y,p x . The covariant vector v 1 , which 
spans the unstable manifold, is indicated in green. Similarly, the covariant vector v 4 
spanning the stable manifold is shown in blue. Of course, these projected vectors are 
not of unit length. The figure may provide an intuitive understanding of the directions 
of maximum stretching or contraction of perturbations in phase space. 

The time evolution of the polar components of v 4 is depicted in Fig. |2j This vector 
was chosen for display since it is the most time consuming to compute and has a non- 
vanishing global exponent. The red lines are the results of a direct integration of Eqs. 
( |2~TT) within the polar environment, whereas the green lines are obtained by converting 
the Cartesian covariant vector Vq to its polar representation. As is indicated in Eq. 
( |24|) . the sign is irrelevant since only the sense of direction counts. The same agreement 
is also obtained for the other covariant vectors (not shown). 

The finite-time Lyapunov exponents (including the limiting case of local exponents) 
in the Cartesian and polar representations are also intimately connected. Let us consider 
a time element of duration t n+1 —t n = r. The Cartesian covariant vector VQ(t n ) evolves 
according to Eqs. (|3]) and ()5]) into a vector exp < A^ r > v^,(t n+ i), which - according 
to Eq. ( )23|) - may be transformed to the polar representation to give 



exp(Ag C0V r)>l(t n+1 )^(Wi: 
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Figure 2. (Color online) Time evolution along the trajectory segment depicted in 
Fig. [T]of the four polar components Sr, 8<f>, 5p r , Sp^ (from top left to bottom right) for 
the covariant vector v p spanning the stable manifold of the spring pendulum. The 
smooth (red) lines are obtained by direct integration with polar coordinates, whereas 
the dashed (green) curves were converted from the Cartesian representation. 



Alternatively, v^(t n ) at t n may be first transformed to the polar representation, 
■M(t n )vc(t n ), and then be evolved in time to the end of the interval, which yields 
the vector 

\M(t n )v l c (t n )\\v £ P (t n+1 ). 



cxp 



Af ov r 



Equating these two expressions and taking the absolute value on both sides yields 



A' 



,cov 



A £,COV + 



In 



\M(t T 



V 



n+1) u C\ i n+l 



(*» 



T 



\M{t n )v l c {t n )\\ ' (25) 

This provides the relation between the covariant local exponents in the two coordinate 
systems. Although expressed here in terms of the Cartesian and polar coordinates, the 
expressions ( |24|) and ( |25l) are completely general. 

Eq. (j25|) clearly shows that the local exponents differ for different coordinate 
systems [3]. To test the last relation, we show by the smooth (red) lines in Fig. [3] the 
local (time-dependent) covariant exponents in the polar representation for the unstable 
manifold (£ = 1, left panel) and for the stable manifold (£ = 4, right panel) of the spring 
pendulum. They are obtained by direct integration of the motion equations (liZT|) with 
polar coordinates. For comparison, the corresponding Cartesian exponents are plotted 
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Figure 3. (Color online) Comparison of time dependent local covariant exponents for 
the unstable (£ = 1, left panel) and stable (£ = 4, right panel) manifolds of the spring 
pendulum. The smooth (red) lines are obtained by numerical integration of the polar 
equations of motion, the dotted (blue) lines by integration of the Cartesian equations 
of motion, followed by a conversion to the polar representation. The agreement is very 
good. For comparison, the corresponding Cartesian exponents, which are the basis of 
the conversion, are shown by the dotted (green) lines. 



by the dashed (green) lines. If the latter are transformed to the polar representation, 
the dotted (blue) lines are obtained. The perfect match with the smooth (red) lines of 
the strictly polar approach asserts the validity of Eq. ( 1251) . 

In Fig. H] we show the (linear) growth rate of the covariance matrix for the 
logarithmic expansion tA^(t) as a function of the averaging interval r for various 
representations. For the GS-vectors this rate matrix is given by [20] 



d2 S (t) 



lim 

T<00 



A? S (r)A? S (r; 



XfXf 



T. 



(26) 



where we still have to distinguish the cases with Cartesian or polar coordinates. For 
covariant vectors the definition is given by Eq. (Q. For clarity we restrict ourselves 
to matrix elements, which do not involve any of the vanishing Lyapunov exponents A2 
and A3. The panel on the right-hand side is a magnification of the small-r regime of 
the left panel. The elements for Gram-Schmidt FTLEs with Cartesian coordinates are 
indicated by red full points (case GS-C), those with polar coordinates by the blue full 
squares (case GS-P). Similarly, the elements involving covariant FTLEs with Cartesian 
coordinates are indicated by green crosses (case COV-C), those with polar coordinates 
by violet open circles (case COV-P). The following observations are made: 

(i) There exists the general symmetry -Dn(r) = D iA {j) and D^ir) = D^iij) for any 
coordinate system and for any vector set used, be it Gram-Schmidt or covariant. 
Wheras the latter equality is trivial, the former is not in view of the fact that neither 
v 1 ,!) nor g x ,g are simply related to each other. 

-Dii(-r) for GS-C and COV-C agree for all r, and similarly D u (r) for GS-P and 
COV-P. This is to be expected, since the covariant vector vi and the GS vector g\ 
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associated with the maximum exponent are identical by construction. 
For small averaging intervals r, the fluctuations of the finite-time exponents as 
measured by Dn(r) and Dn{r) significantly differ for the cases GS-C and GS-P, as 
has been observed before [22]. A similar observation is made for the cases COV-C 
and COV-P. These differences gradually disappear for averaging times r > 10. 
The cross correlation D\±(t) for the covariant and Gram-Schmidt vector sets differ 
enormously for small r, regardless what coordinate system is used. The difference 
is even larger between COV-P and GS-P than between COV-C and GS-C. By 
increasing r these differences very slowly disappear. 

From the left panel of Fig. H] one infers that for r > 200 all elements of the 
covariance matrix for the FTLEs become independent of the choice of the set of 
perturbation vectors and of the coordinate system. Thus, lim^oo D{r) becomes 
a dynamical invariant independent of the metric and parameterization |20j. For 
high- dimensional chaotic systems such a result has been interpreted as evidence 
for effective hyperbolicity of the dynamics and of the statistical insignificance of 
hyperbolic tangencies. 
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Figure 4. (Color online) Chaotic pendulum: r dependence of the growth rate of the 
covariance matrix for the logarithmic expansion tA^(t). Only matrix elements not 
involving vanishing global exponents are shown. The panel on the right-hand side is 
a magnification of the small r regime. Red full points: Cartesian coordinates and GS 
exponents (case GS-C); blue full squares: polar coordinates and GS exponents (case 
GS-P); green crosses: Cartesian coordinates and covariant exponents (case COV-C): 
violet open circles: polar coordinates and covariant exponents (case COV-P). In all 
cases the evolution is in the direction of the positive time axis. 



4. The Henon-Heiles system 



As a second example, we consider the familiar Henon-Heiles system [21] with a 
Hamiltonian in Cartesian coordinates 

lipl + pI) + l(x 2 + y 2 ) + % 2 y - lv 3 - (27) 



h, 



c 
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For an energy 1/6, the system is known to be chaotic (with a global Lyapunov spectrum 
{0.127 7 , 0, 0, — 0.127 7 }), where the trajectory visits most of the accessible phase space 
[251 I2T] . The local Lyapunov exponents display the symplectic pairing symmetry Eq. 
(IIOHIII) . and the time reversal symmetry of Eq. ( 1T3|) . as is demonstrated in the respective 
left and right panels of Fig. |5j 
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Figure 5. (Color online) Local Lyapunov exponents for the Henon-Heiles system with 
energy H = 1/6. Left panel: Illustration of the symplectic local pairing symmetry, 
Eq. (|10lll|) . for the Gram-Schmidt (GS) exponents ( + ^Ap S and (+)Ap S (smooth red 
lines). The dimension of phase space D — 4. The inequality Eq. (fT2)) is demonstrated 
by the dashed green line for (~)Ag . Right panel: Verification of the time-reversal 
invariance property (1131) for the local covariant exponents specified. 



In polar coordinates defined by x = r cos <fi and y = r sin <fi the Hamiltonian becomes 

(28) 



tt Vl , V\ r- r- , 

tip = ^r H sin 3<p, 

F 2 2r 2 2 3 



where, as before, p r = f and p^ = r 2 <p are the conjugate momenta. With an analysis 
completely analogous to that in the previous section for the spring pendulum, the 
covariant vectors in the polar representation may be obtained. Here, it suffices to 
present in Fig. [H] only results for the polar covariant vector Vp associated with the 
maximum exponent. One observes that the covariant vectors are unique (up to the 
sign as the example in the previous section testifies), and that it does not matter what 
parameterization is used for their computation. 

5. Resume 



There are many recent applications of Lyapunov vectors, which range from weather 
forecasting [28J, geophysical applications [27], and studies of transport in turbulent flow 
[29], [30]. to identifying and probing rare chaotic events by Lyapunov- weighted dynamics 
[31] and/or path sampling [32J. 

From the two sets of perturbation vectors commonly used - the orthonormal Gram- 
Schmidt vectors and the co-moving covariant vectors - the latter are more attractive 
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Figure 6. (Color online) Time evolution of the four polar components Sr, 6<j>, 5p r , Sp^ 
(from top left to bottom right) for the covariant vector vp spanning the unstable 
manifold of the Henon-Heiles system. In every panel two lines are shown. The smooth 
(red) lines are obtained by direct integration with polar coordinates, whereas the 
dashed (green) curves were converted from the Cartesian representation. Both lines 
agree and cannot be distinguished. 



from a physics point of view. They reflect time-reversal invariance and, as is shown 
here, they are unique and independent of the choice of the coordinate system. By 
uniqueness we mean that they are properties of the flow in tangent space, and when 
they are known in a particular frame, they may be easily converted to another. They 
provide a spanning set for the Oseledec splitting of the tangent space into a hierarchy 
of subspaces, each characterized by a well-defined Lyapunov exponent. In this sense, 
they are a property of the physical system. This is not the case for the more artificial 
ortho-normal Gram-Schmidt vectors, which are a mathematical device for computing 
volume changes of d- dimensional volume elements, d < D, in the .D-dimensional phase 
space. Their advantage is that they are easier and faster to compute. According to 
the algorithm of Ginelli et al. [TT], the computation of the covariant vectors at a 
phase-point T requires the knowledge of the Gram-Schmidt vectors for all points on the 
trajectory through T both in the past and the future j2]. Due to the availability of fast 
computers and programming techniques, rather high- dimensional systems have already 
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been examined jU |33J EH] 

For many problems only the first Lyapunov vector associated with the maximum 
(global) exponent is required. Since the first Gram-Schmidt and the first covariant 
vectors always agree, one gets away with the fast computation of the first Gram-Schmidt 
vector, which is covariant. This advantage may dissolve again if the application requires 
the knowledge of the maximum local exponent together with its covariant vector. Due 
to possible exponent entanglement [T3| 116] the maximum local exponent at a particular 
instant may belong, say, to the fifth vector. In such a case more than one covariant 
vectors need to be computed. 
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